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Nonequilibrium spin texture within a thin layer below the surface of current-carrying 
topological insulator 612803: A first-principles quantum transport study 
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We predict that unpolarized charge current injected into a ballistic thin film of prototypical 
topological insulator (TI) Bi2Se3 will generate a noneollinear spin texture S(r) on its surface. Fur¬ 
thermore, the nonequilibrium spin texture will extend into ~ 2 nm thick layer below the TI surfaces 
due to penetration of evanescent wavefunctions from the metallic surfaces into the bulk of TI. Av¬ 
eraging S(r) over few A along the longitudinal direction defined by the current flow reveals large 
component pointing in the transverse direction. In addition, we find an order of magnitude smaller 
out-of-plane component when the direction of injected current with respect to Bi and Se atoms 
probes the largest hexagonal warping of the Dirac-cone dispersion on TI surface. Our analysis 
is based on an extension of the nonequilibrium Green functions combined with density functional 
theory (NEGF+DFT) to situations involving noneollinear spins and spin-orbit coupling. We also 
demonstrate how DFT calculations with properly optimized local orbital basis set can precisely 
match putatively more accurate calculations with plane-wave basis set for the supercell of 612803. 

PACS numbers: 72.25.Dc, 75.70.Tj, 71.15.Mb, 72.10.Bg 


The newly discovered three-dimensional topological in¬ 
sulator (3D TIs) materials possess a usual band gap in 
the bulk while also hosting metallic surfaces. The low- 
energy quasiparticles on these surfaces behave as mass¬ 
less Dirac fermions whose spins are locked to their mo¬ 
menta due to strong spin-orbit coupling (SOC).® Such 
spin-momentum locking is viewed as a resource for spin- 
tronic applications.!^ For example, very recent experi¬ 
ment have demonstrated magnetization dynamics of a 
single ferromagnetic metallic (FM) overlayer deposited 
on the surface of 3D TIs due to current-induced SO 
torques. Another recent experimenfPI has detected spin- 
to-charge conversiorP^ when precessing magnetization of 
the FM overlayer pumps pure spin current into the metal¬ 
lic surface of 3D TIs. 

The microscopic mechanism behind these phenomena 
can be traced to the so-called Edelstein effect (EE), orig¬ 
inally predictecP for a diffusive two-dimensional elec¬ 
tron gas (2DEG) with the Rashba SOCP and observed 
much later experimentally.!^ In the EE in 2DEG, longi¬ 
tudinal unpolarized charge current flowing along the x- 
axis drives a homogeneous nonequilibrium spin density 
S = (0,5'y,0) pointing in the transverse direction. The 
diffusive metallic surface of TIs also exhibits EE where a 
current-driven spin density S is substantially enhancecP 
(by a factor hop jan ^ 1, with vp being the Fermi veloc¬ 
ity in TI and Oip is the strengtiP of the Rashba SOC in 
2DEG). This enhancement can be explained by the spin- 
momentum locking along the single Fermi circle,!!^ formed 
in /c-space at the intersection of the Dirac cone energy- 
momentum dispersion and the Fermi energy plane, in 
contrast to spin-momentum locking along the two cir¬ 
cled in the case of Rashba 2DEG which counter the 
effect of each other. This has motivated recent exper¬ 
iment^^ probing S directly in three-terminal geometry 
where nonmagnetic electrodes inject unpolarized charge 
current into a TI, while a third FM contact deposited 



FIG. 1. Schematic view of a two-terminal setup where a thin 
film of Bi2Se3 is attached to two macroscopic reservoirs biased 
by the electrochemical potential difference eVh — [il — /xr. 
The elean Bi2Se3 film is infinite along the a:-axis (i.e., the di¬ 
rection of transport) and the y-axis, while its thickness along 
the z-axis is chosen as 5 QLs. The shaded cell of length 
dx — 5 A defines volume for averaging S(r) from Fig.[^ which 
is then plotted in Fig. over the corresponding cross section 
of thin film within the yz-plane. 


in the middle of the top surface of the TI film detects a 
voltage signal when a non-zero S is induced. These se¬ 
tups quantify the projection of S onto the magnetization 
of the third FM contact. 

However, this picture of EE on the surface of TI is 
based on simplistic model Hamiltonians.!^ Here we em¬ 
ploy first-principles quantum transport approach to an¬ 
alyze microscopic details, over ^ 1 A length scale, of 
current-driven S(r) in the two-terminal ballistic thin film 
geometry hosting realistic TI material, as illustrated 
in Fig. [2 We choose Bi 2 Se 3 as the prototypical TI 
material—with its single Dirac cone in the surface band 
structure, relatively large bulk band gap, and Dirac point 
(DP) inside the gap [see Fig. |^a)]— on which many re¬ 
cent experiments probing EE directl}ff^ or indirectl}®! 
have been ^rformed. The central region of the de¬ 
vice in Fig. [ij which has length Lx = 21.5 A along the 
x-axis and infinite width along the ^-axis, is attached 
to two semi-infinite electrodes made of the same mate- 
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FIG. 2. (a) The arrangement of Bi and Se atoms within a supercell of Bi 2 Se 3 thin film of thickness 5 QLs. The inset in panel 
(a) shows BZ in the kx-ky plane at = 0 with special /c-points T, M, and K marked, (b) The vector field of nonequilibrium 
S(r) within selected planes shown in (a), generated by injection of unpolarized charge current along the a:-axis (see also Fig.[^. 
The planes 1 and 3 correspond to the top and bottom metallic surfaces of Bi 2 Se 3 thin film, while plane 2 resides in the bulk at 
a distance d ~ 0.164 nm away from plane 1. (c) The vector fields in (b) projected onto each of the selected planes in (a). The 
real space grid of r-points in panels (b) and (c) has spacing ~ 0.4 A. 


rial. The electrodes are assumed to terminate at infin¬ 
ity into macroscopic Fermi liquid reservoirs where elec¬ 
trons are thermalized to acquire electrochemical poten¬ 
tial jiiL in the left reservoir and jur in the right one. The 
Hamiltonian of the central region and the electrodes is 
obtained from the noncollinear density functional the¬ 
ory (ncDFT), implemented by us in ATK packag^^l^ 
ing optimized pseudo-atomic localized basis functionals 
and SOC introduced via the total-angular-momentum- 
dependent pseudopotentials. IIS The transport properties 
of the system in Fig. are computed using the nonequi¬ 
librium Green function (NEGF) formalism, IIS so that 
our approach represents an extension of the widely used 
NEGF-}-DFT frameworlilS to transport problems involv¬ 
ing noncollinear spins and SOC. 

In the simplest picture—based on effective Hamilto¬ 
nian = vf{^ X p) ■ Gz (S' is the vector of the Pauli 
matrices; p is the momentum operator; and Sz is the 
unit vector along the z-axis in Fig. describing mass¬ 
less Dirac electrons on the metallic surfaces of TIs—the 
spin and momentum of electronic eigenstates are orthog¬ 
onal to each other along the single Fermi circle. This 
generates net homogeneous S = (0,5'y,0) after a n ap¬ 
plied electric field shifts the Fermi circl ej^ l ^H ^ along 
the momentum parallel to E^. Such man ifestation of 
EE persists in ballistic samples as welPm^ where there 


is no electric field within the TI but instead one applies 
bias voltage eVi, = /i^ — to inject a current into the 
TI, as illustrated in Eig The relations Sy oc Ex or 
Sy OC Vb describing EE in the diffusive or ballistic trans¬ 
port regimes, respectively, are allowed only in nonequi¬ 
librium since in equilibrium S changes sign under time 
reversal, and, therefore, has to vanish (assuming absence 
of magnetic field). 

This simplistic picture can be contrasted with our prin¬ 
cipal results in Eigs. and When a small (ensur¬ 
ing linear-response transport regime) V 5 is applied be¬ 
tween the reservoirs in Eig.[^ the unpolarized charge cur¬ 
rent injected into Bi 2 Se 3 thin film generates a nonequi¬ 
librium S(r) whose complex noncollinear texture within 
three planes selected in Eig. [^a) is plotted in Eigs. |^b) 
and [^c). Eor the visualization we use real-space grid 
for r whose spacing is 0.4 A. Eurthermore, Eigs. 
and [^demonstrate that nonequilibrium spin texture S(r) 
will appear not only on the TI surface, but also within 
^ 2 nm thick layer of its bulk just below the top and 
bottom surfaces. This feature is explained in Eigs. ia) 
and [^b) showing spatial profile of the local density of 
states (LDOS) at the Eermi energy Ep over the cell de¬ 
picted in Eig. The non-zero LDOS and the corre¬ 
sponding S(r) in the bulk of the TI thin film stem from 
evanescent wavefunctions which originate from the top 
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and bottom metallic surfaces and penetrate into the en¬ 
ergy gap of the insulating bulk. The Bi 2 Se 3 is a strongly 
anisotropic material composed of quintuple layers (QLs) 
of Bi and Se atoms, illustrated in Fig. [^a), where one 
QL consists of three Se layers strongly bonded to two Bi 
layers in between. For Bi 2 Se 3 film thinner than 5 QLs, 
the evanescent wavefunctions from the top and botto m 
metallic surface can overlap to create a minigapP^^ at 
the DP. We select the thickness of Bi 2 Se 3 to be 5 QLs 
along the 2 ;-axis in Fig. which ensures that the LDOS 
in Figs, [^a) and [^b) goes to zero on the plane half way 
between the top and bottom surfaces of the TI thin film. 

Upon averaging nonequilibrium S(r) over a 5 A 
long cell depicted in Fig. we obtain spatial profiles in 
Figs. j^c) and id) which show that Sy is the largest 
component independently of the direction of incoming 
electrons. An order of magnitude smaller Sz component 
shown in Fig. [^d) appears for electrons incoming along 
current direction 2 marked in panel (e). This is in ac¬ 
cord with experiments in equilibrium where spin- and 
angle-resolved photoemission spectroscop}!^ finds largest 
out-of-plane spin component along the corresponding di¬ 
rection in the 2D Brillouin zone (BZ). This is due to 
hexagonal warping of the Dirac cone s urface band struc¬ 
ture, as confirmed by DFT calculation^^^^^ finding that 
equilibrium expectation value of spin in the eigenstates 
of Bi 2 Se 3 surfaces tilts out of the 2D BZ. Thus, Fig. 
offers a novel prescription for probing hexagonal warp¬ 
ing even close to DP via transport measurements where 
charge current is injected in different directions relative 
to the orientation of the lattice of Bi and Se atoms. 

We now explain the technical details of our calcu¬ 
lations. The extension of DFT to the case of spin- 
polarized systems is formally derived in terms of total 
electron density n(r) and vector magnetization density 
m(r). In the collinear DFT, m(r) points in the same 
direction at all points in space, which is insufficient to 
study magnetic systems where the direction of the local 
magnetization is not constrained to a particular axis or 
systems governed by SOC. In ncDFT,^ the exchange- 
correlation (XC) functional £’xc['^(i *)5 depends on 

m(r) pointing in arbitrary directions. The local den¬ 
sity approximation(LDA) and most often employed ver¬ 
sion of generalized gradient approximation (GGA), im¬ 
plemented also by us in ATK^ make additional ap¬ 
proximation^^ that lead to the XG magnetic field 
Bxc(r) = (5F^xc[^(r), m(r)]/(im(r) being parallel every¬ 
where to m(r). 

The single-particle spin-dependent Kohn-Sham (KS) 
Hamiltonian in ncDFT takes the form 

^KS = + ^H(r) + Uxc(r) + Uext(r) - tr • Bxc(r), 

( 1 ) 

where Uext(i*) and Uxc(i*) = 

F^xc[^(i*), m(r)]/(in(r) are the Hartree, external and XG 
potential, respectively. Diagonalization of ^ks proceeds 
by approximating the Hilbert space of all single-electron 
eigenfunctions with a finite set of basis functions. A 
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FIG. 3. (a),(b) The spatial profile of LDOS at Ep over the 
cross section of the cell denoted in Fig.[^ (c),(d) The spatial 
profile of the components of S(r), obtained by averaging its 
texture plotted in Fig. over the cell of Bi 2 Se 3 thin film 
marked in Fig.^ The direction of injected charge current for 
the results in panels (a) and (c), or the results in panels (b) 
and (d), is denoted in panel (e) relative to the orientation of 
the lattice of Bi and Se atoms. The bottom surface of Bi 2 Se 3 
is located at z = 0 nm, and the top TI surface is located at 
4.56 nm. 


popular basis set is plane-waves (PWs), where varying 
only one parameter (the energy cutoff) allows one to 
improve the basis systematically. Linear combination of 
atomic orbitals (LGAO) basis sets require more tuning, 
however, they simplify the NEGF calculation^^ where 
one has to spatially separate system into the central 
region and semi-infinite electrodes, as illustrated in 
Fig. [2 

Since the pioneering screeninj^ of candidate TI ma- 
terials via ncDFT calculations, their elec tronic band 
structure has most often been calculatecP^l^ using 
PW ncDFT with electron-core interactions described 
via projector augmented wave (PAW) method.^ In 
Fig. I^a) we demonstrate that such calculations, per- 
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FIG. 4. (a) The electronic band structure computed for a 

supercell of Bi 2 Se 3 film shown in Fig. [^a) using LCACp^ 
ncDFT implemented in ATK?^ package. This is compared 
with the electronic band structure obtained using PW ncDFT 
implemented in VASP package.^ (b) Zero-bias transmission 
function of Bi 2 Se 3 thin film in the two-terminal geometry of 
Fig. for electrons injected along the F-M direction {ky = 
0) in the inset of Fig. [^a), computed using NEGF+ncDFT 
formalism implemented in ATK package. 

formed by VASP package}^ can be accurately repro¬ 
duced by pseudopoteutial-based LCAO ncDFT imple¬ 
mented in ATK.^ The super cell considered in both cal¬ 
culations is shown in Fig. |^a), which includes 5 QLs 
terminated by Se atomic layer on both the top and bot¬ 
tom surface, as well as 7.5 A thick vacuum layer above 
and below these Se atomic layers. 

We note that previous attempt^ to apply 
pseudopotential-based LCAO ncDFT to Bi 2 Se 3 have 
yielded either poor accuracy of its electronic band 
structure (e.g., compare our Fig. [^a) with Fig. 1 
in Ref. [2T]) or have required intricate fine tuning 
Therefore, we provide here a complete recipe for the 
proper usage of LCAO ncDFT to reproduce Fig. |^a). 
In ATK calculations in Fig. [^a), the electron-core 
interactions are described by norm-conserving pseudopo¬ 
tentials. The pseudopotentials are obtained by mapping 
the solution of the Dirac equation, which naturally 
includes SOC,^!^ to non-relativistic pseudopotential, 
^ps = Vl + “b with local contribution W 

and non-local contributions Vnl from the total angular 
momentum j = 1 + 1/2 and j = I — 1/2. The non-local 
terms are expanded in terms of SO projector functions, 

++ = , where *^(± 1 / 2 ,^ are 

normalization constants and the indices (a,/3 denote the 
possible spin orientations We use Perdew-Burke- 

Ernzerhof (PBE) parametrization of GGA for the XC 
functional and a LC AO ba sis set {0^} generated by 
the OpenMX package ,1^^^ which consists of s2p2dl 


orbitals on Se atoms and s2p2d2 on Bi atoms. These 
pseudoatomic orbitals were generated by a confinement 
schem^I^ with the cutoff radius 7.0 a.u. and 8.0 a.u. 
for Se and Bi atoms, respectively. The energy mesh 
cutoff for the real-space grid is chosen as 75.0 Hartree. 
In VASP calculation^^ in Fig. |^a), the electron-core 
interactions are described by PAW method,!^ and we 
employ PBE GGA for the XC functional. The cutoff 
energy for the PW basis set is 350 eV. In both ATK and 
VASP calculations we employ 11 x 11 x 1 /c-point mesh 
within Monkhorst-Pack scheme for the BZ integration. 

The eigenstates of the KS Hamiltonian in 

Eq. 0 make it possible to construct the equilibrium 
density matrix Peq = electrons 

at jiL = I^R and temperature T described by the 
Fermi distribution function f{E). The local electron 
and magnetization density, as the central variables of 
ncDFT, are obtained from n(r) = (r|Trspin[Peq] A) and 
m(r) = (r|Trspin[peqCr] |r), where the trace is taken over 
the spin Hilbert space. 

In steady-state nonequilibrium due to dc current flow¬ 
ing between the left and right reservoirs in Fig. we 
construct the nonequilibrium density matrid^ pneq us¬ 
ing NEGFs: 

-|-oo 

Pneq = J dE (E) — Peq. (2) 

— oo 

This yields S(r) = |(r|Trspin[PneqO’]|r) plotted in 

Figs. and The NEGF formalisnP^l for steady- 

state transport operates with two central quantities— 
the retarded GF, G{E), and the lesser GF, G^{E )— 
which describe the density of available quantum 
states and how electrons occupy those states, re¬ 
spectively. In the absence of inelastic processes, 

these are given by G = [EO — Hks — ^l — and 

G< = iG[fLTL + /i?r_R]G^ Here the self-energies 'El,r 
are due to semi-infinite electrodes, fh.R = f{E — I^l^r) 
and = i{EtL,R ~ level broadening ma¬ 

trix. For the chosen LCAO basis set, the Hamilto¬ 
nian matrix Hks is composed of elements (0i|^Ks|0j) 
and the overlap matrix O is composed of elements 
{(j)i\(j)j). In the linear-response transport regime con¬ 
sidered here, Eq. § can be expandecP^ to first or¬ 
der in bias voltage 14. Since S(r) is zero in equilib¬ 
rium (because of assumed absence of external magnetic 
field), the linear-response density matrix can be sim- 

-|-CXD , 

plified toPlpneq = ^ / dEGTiGt (-|^). Other- 

wise, the gauge-invariant form of pneq requires additional 
term^^ to properly remove the equilibrium expectation 
value of a considered physical quantity. 

The retarded GF also allows us to obtain the 
transmission function of the device in Fig. 
T{E, ky) = Tr[ri^GrLG’*‘], which depends on energy and 
transverse momentum ky due to assumed periodicity in 
the ^-direction. The total transmission function T{E) 
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is obtained by integrating over ky^ which determines the 
linear-response conductance via the Landauer formula, 

G = ^JdET{E) (-§£)■ We confirm in Fig.j^b) that 

T{E^ky = 0) = 2 for ^ within the bulk gap shown 
in Fig. [^a) because only one doubly degenerate helical 
conducting channel is open for transport in that energy 
rang^^ for injected electrons with momentum along the 
F-M direction {ky = 0). 

In conclusion, using NEGF+ncDFT framework imple¬ 
mented by us in ATK package,^ we computed a nonequi¬ 
librium spin texture S(r) within a thin film of current- 
carrying Bi 2 Se 3 TI material. The non-zero texture ap¬ 
pears on the TI metallic top and bottom surfaces, as 
well as within bulk layers of thickness 2.0 nm below 
the surfaces that effectively dope the bulk by evanescent 
wavefunctions. The spin texture is noncollinear and com¬ 
plex on length scales < 1 A. Upon averaging it over a few 
A we find a simpler pattern—with either S = (0,5'y,0), 


or S = where Sy/S'^ ^ 1—depending on the 

direction of injected current with respect to orientation 
of the lattice of Bi and Se atoms. Such dependency offers 
a novel probe, via electronic transport measurement s,l^ 
of the hexagonal warping of the Dirac cone surface band 
structure. For the envisaged spintronic applications of 
TIs, it is essential to understand how S(r) changes due to 
finite bias voltage or self-consistent couplin^^ to magne¬ 
tization of a ferromagnetic (metal or insulator) overlayer, 
which we relegate to future studies. 
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